Hepatic transcriptome analysis identifies genes, polymorphisms and pathways involved in the fatty acids metabolism in sheep

Fatty acids (FA) in ruminants, especially unsaturated FA (USFA) have important impact in meat quality, nutritional value, and flavour quality of meat, and on consumer’s health. Identification of the genetic factors controlling the FA composition and metabolism is pivotal to select sheep that produce higher USFA and lower saturated (SFA) for the benefit of sheep industry and consumers. Therefore, this study was aimed to investigate the transcriptome profiling in the liver tissues collected from sheep with divergent USFA content in longissimus muscle using RNA deep-sequencing. From sheep (n = 100) population, liver tissues with higher (n = 3) and lower (n = 3) USFA content were analysed using Illumina HiSeq 2500. The total number of reads produced for each liver sample were ranged from 21.28 to 28.51 million with a median of 23.90 million. Approximately, 198 genes were differentially regulated with significance level of p-adjusted value <0.05. Among them, 100 genes were up-regulated, and 98 were down-regulated (p<0.01, FC>1.5) in the higher USFA group. A large proportion of key genes involved in FA biosynthesis, adipogenesis, fat deposition, and lipid metabolism were identified, such as APOA5, SLC25A30, GFPT1, LEPR, TGFBR2, FABP7, GSTCD, and CYP17A. Pathway analysis revealed that glycosaminoglycan biosynthesis- keratan sulfate, adipokine signaling, galactose metabolism, endocrine and other factors-regulating calcium metabolism, mineral metabolism, and PPAR signaling pathway were playing important regulatory roles in FA metabolism. Importantly, polymorphism and association analyses showed that mutation in APOA5, CFHR5, TGFBR2 and LEPR genes could be potential markers for the FA composition in sheep. These polymorphisms and transcriptome networks controlling the FA variation could be used as genetic markers for FA composition-related traits improvement. However, functional validation is required to confirm the effect of these SNPs in other sheep population in order to incorporate them in the sheep breeding program.


Introduction
Meat quality is an economically important trait because of consumer's choice which includes both visual and sensory traits, health benefits, and humane production system. Recently, fatty acids (FA) composition is being considering as a new feature for lamb quality [1]. Ruminants' meat is generally containing higher levels of saturated fatty acids (SFA), which are widely correlated with health problem such as heart disease, stroke, and obesity [2], so consumers are favouring leaner meats containing less SFA and higher polyunsaturated fatty acids (PUSFA) [3,4]. PUSFA, mainly omega-3 are considered beneficial for human health that reduce the serum low density lipoprotein (LDL)-cholesterol, total cholesterol concentration, and modulate immune functions [5]. Additionally, desirable sensorial characterisctic of meat is associated with PUSFA and MUSFA (monounsaturated fatty acids) [6]. Note, sheep meat is rich in omega-3 long-chain (�20) FA (ω3 LC-PUSFA), eicosapentaenoic (EPA, 20:5ω3), and docosahexaenoic (DHA, 22:6ω3) which are beneficial for human health and immunity [7]. Meat production with a higher PUSFA and lower SFA content is, therefore, important to improve human health without requiring substatial changes in customers' habit of meat consumption.
Molecular breeding is recommended as one of the most realistic approaches for increasing PUSFA-and reducing SFA-content. However, identification of the candidate genes and genomic networks is the first step to achieve the goal. Notably, FA compositions are the welldefined compounds describing the phenotypic traits which are possible to improve through genetic selection. FA compositions show moderate to high heritability ranging from 0.15 to 0.63 [8,9]. Identification of genetic factors controlling FA composition could be implemented in breeding programmes to select animals that produce higher PUSFA and lower SFA in meat. Therefore, it is crucial to understand the genomics of FA metabolism to select sheep with higher PUSFA and lower SFA content. FA metabolism is a complex process, which involves lipolysis of dietary fat, biohydrogenation in the rumen, and de novo synthesis of FA by rumen bacteria. Furthermore, absorption and transport of FA by the host animal, de novo synthesis, elongation and desaturation in the animal's tissues, hydrolysis of triglycerides, esterification, and the oxidation of FA or its metabolization into other components together make it a complex process to decipher [10].
High-throughput sequencing technologies (RNA-Seq) are now widely using for transcriptome analysis because of an unprecedented accuracy and data insight [11]. The reliable and comprehensive data from RNA-Seq can not only describe the genes' structure, but also provide a better understanding of the biological function of genes [12]. This technology is allowing the animal breeding industry to significantly increase the rate of genetic progress [13]. Several recent studies have used RNA deep sequencing to identify differentially expressed genes related to FA metabolism in muscle and liver in domestic animals such as in pigs [14,15], and cattle [16]. But our understanding of genomic signature behind the FA metabolism in sheep at the molecular level is limited. Although several candidate genes, such as ACACA [17], FASN and SCD [18] are reported to be associated with FA and fat content in various sheep breeds, the whole genomics underlying the FA metabolism in sheep is remained to be deciphered. In accordance with other studies of FA composition, there is an inevitable need for using RNA deep sequencing for transcriptome profiling related to higher PUSFA and lower SFA in sheep. Therefore, the aim of this study was to elucidate the genes and pathways involved in FA metabolism in the liver tissue using RNA deep sequencing technology. For this purpose, differential expression analysis of transcriptome was performed in the liver tissues collected from sheep with higher and lower USFA in their longissimus muscle. In addition, gene polymorphism and association analyses were also performed for the putative candidate genes. Since consumers intake FA from muscle tissues, the longissimus dorsi muscle tissues were used for FA composition analysis; whereas FA are metabolised in the liver so hepatic transcriptome analysis was performed to unravel the genes and networks controlling FA metabolism in sheep.

Quality control and analysis of RNA deep sequencing data
From the sheep (n = 100) population, liver tissues with higher (n = 3) and lower (n = 3) unsaturated fatty acids (USFA) content were selected for high-throughput sequencing. cDNA libraries from 6 samples of sheep liver tissues (3 from HUSFA = higher USFA, and 3 from LUSFA = lower USFA) were sequenced using Illumina HiSeq 2500. The sequencing produced clusters of sequence reads with maximum of 100 base-pair (bp). After quality control and filtering, the total number of reads for liver samples were ranged from 21.28 to 28.51 million with a median of 23.90 million. Total number of reads for each group of samples and the number of reads mapped to reference sequences are shown in Table 2. In case of LUSFA group, 84.51 to 85.69% of total reads were aligned to the reference sequence, whereas 85.20 to 87.38% of the total reads were aligned in case of the HUSFA group.

Differential gene expression analysis
Differential gene expression from livers tissues of sheep with HUSFA and LUSFA levels were calculated from the raw reads using the R package DESeq. The significance scores were corrected for multiple testing using Benjamini-Hochberg correction. A negative binomial distribution-based method implemented in DESeq was used to identify differentially expressed genes (DEGs) in the liver tissues collected from sheep with divergent unsaturated fatty acids (USFA) level in the longissimus muscle. A total of 198 DEGs were selected from the differential expression analysis using criteria p adjusted < 0.05 and log2 fold change > 1.5 (Fig 1). In liver tissues, 110 genes were found to be highly expressed in HUSFA group, whereas 98 genes were found to be highly expressed in LUSFA group (S1 Table). The range of log2 fold change values for DEGs were between 4.09 to-4.80 (Fig 2 and Table 3). Heatmaps illustrated the top 30 upand down-regulated genes identified in the liver tissues from HUSFA and LUSFA sheep. The top 30 up-and down-regulated genes identified in the liver tissues with divergent USFA levels along with log FC and p values are listed in the Table 3. The differential expression analysis of data revealed both novel transcripts and common genes which were previously identified in various gene expression studies related to FA. Novel transcripts from this analysis and commonly found genes are mentioned in detailed in the discussion section.

Biological function analysis for DEGs
Gene ontology (GO) and pathway enrichment analysis were performed to gain insight into the predicted genes networks. The most significant GO terms were categorized into biological processes, cellular components, and molecular functions (Fig 3). The enriched biological processes identified were mainly related to cytokinesis, glycoprotein metabolic process, mitotic spindle, N linked glycosylation, acute inflammatory response, and regulation of developmental

PLOS ONE
process. Cellular components consisted of cell projection part, extracellular space, integral to plasma membrane, and proteinaceous extracellular matrix were significantly enriched. The molecular functions identified were related to kinase inhibitor activity, growth factor binding, and GTPase activity. A total of 11 significantly enriched KEGG pathways were identified as overrepresented for the DEGs. The KEGG pathway analysis showed that glycosaminoglycan biosynthesis-keratan sulphate, adipokine signaling, galactose metabolism, endocrine and other factor-regulated calcium metabolism, mineral metabolism, and PPAR signaling pathways were significantly involved in fatty acids metabolism regulation in the liver (Fig 4).

Regulatory hub genes of the hepatic transcriptome network
In order to identify the key regulatory genes in the transcriptional network, a liver-specific protein-protein interaction (PPI) network was created that comprised of 48 seed genes; and 530 nodes connected with 578 edges. Based on the network centrality measures, the potential Hub genes were identified, among which SOCS3, CBX6, MCM4, ITGB3, TGFBR2, GPRASP1, CELSR3, SDC3, SPOCK1, SEL1L and LEPR were upregulated, whereas ACTA2, GPRASP1, TPM2, TGM3, PTK6, and LTF were downregulated (Fig 5A and 5B). In addition, we have also  created a liver-specific gene co-expression network to pick up more potential Hub genes, those could have been missed in the PPI network. The co-expression network illustrated that RAC-GAP1, MCM4, SDC3, CKAP2, RNASE6, PREX1, QSOX1, and FUT11 were the upregulated, whereas CDC42EP5, SSC5D, GPRASP1, HRC, NRN1 and TPM2 were the downregulated Hub genes (Fig 6A and 6B). Notably, RACGAP1, TGFBR2, LEPR, MCM4, SDC3, GPRASP1 were the common Hub genes in both PPI and co-expression network analysis (S2 and S3 Tables).

Validation of selected DEGs using quantitative Real Time PCR (qRT-PCR)
A total of 8 differentially expressed genes (CYP17A1, FABP7, GSTCD, SLC25A30, APOA5, GFPT1, LEPR and TGFBR2) were selected and quantified using qRT-PCR, as part of RNA-Seq results validation. For this purpose, the same samples used in the RNA-deep sequencing were used. Comparison of qRT-PCR data for 8 selected genes showed quantitative concordance of expression with the RNA-Seq results (Fig 7). Gene expression values for qRT-PCR were normalized using the average expression values of housekeeping gene GAPDH and β-Actin. Details of GenBank accession numbers, primers sequences, product size, and annealing temperature for qRT-PCR validation used in this study are listed in Table 4.

Gene variation analysis and association study
A total of 226 single nucleotide polymorphisms (SNPs) were identified in 31 DEGs between higher and lower USFA groups (S4 Table). The selected polymorphisms identified in DEGs for liver samples are given in Table 5. The distribution of the number of genes having SNPs, and selected SNPs used for validation are shown in Fig 8A and 8B, respectively. Validation of the SNP results for the association study was carried out by selecting a total of 4 SNPs based on the functional SNPs and the function related to fatty acid metabolism (Fig 8B and S5 Table).
The selected SNPs were harboured in APOA5, CFHR5, TGFBR2 and LEPR genes. These SNPs

PLOS ONE
were analysed to validate their segregation and association in the studied sheep population (n = 100). Our association analyses suggested that, the polymorphisms in APOA5, CFHR5, TGFBR2 and LEPR were associated with fatty acid composition (Table 6) in the studied sheep population.

Analysis of RNA seq data
This study describes the transcriptome profiles of the liver tissues collected from sheep with higher and lower unsaturated fatty acids (HUSFA vs LUSFA) content in their longisimuss muscle. RNA-Seq have allowed for the large-scale analysis of genomic data, providing new opportunities for the characterization of transcriptome architectures [19]. According to the mapping results, the average number of reads was 23.90 million reads, and on an average 85.89% of the reads were categorized as mapped reads corresponding to exon reads ( Table 2).
High-quality reads of mapping results were obtained from an RNA-Seq analysis of the six libraries by comparing to the Ovis aries genome. The proportion of reads mapped to exons of annotated genes was in accordance (85.70-86.95%) with the previous studies [20][21][22] in sheep muscle transcriptome, but was higher than that reported by Wang et al. [12] (68.97%) in short-tailed sheep adipose tissue. The percentage of annotated reads varies from 66.40% to 86.95% in sheep transcriptome studies [12,[20][21][22] supporting our results. The differences between mapping percentages might be due to the current reference transcriptome assembly that might not cover all the transcribed mRNA [23] and consequently low abundant transcripts are less likely to be mapped to the transcriptome assembly [24]. Illumina sequencing data have been described as replicable with relatively little technical variation [25]. Therefore,

Differential express gene analysis
A total of 198 genes were differentially regulated in liver tissues from sheep with divergent USFA levels (S1 Table). The top up-and down-regulated genes in the liver tissues were Zinc Finger Protein 549 with log2 fold change 4.09, and olfactory receptor-like protein DTMT with log2 fold change -4.80, respectively ( Table 3). The genes encode Zinc-finger proteins are involved in cell proliferation and differentiation [26] as well as regulate lipid metabolism [27]. However, the relation between olfactory receptor family genes and USFA is yet to understand. Among the DEGs screened with stringent criteria in the present study, a large proportion of key genes involved in FA biosynthesis, fat deposition, adipogenesis, and lipid metabolism were identified, such as APOA5, SLC25A30, GFPT1, LEPR, TGFBR2, FABP7, GSTCD and CYP17A. APOA5 regulates the assembly and secretion of lipoproteins [28] and controls the plasma triglyceride levels in humans and mice [29,30]. Interestingly four members of SLC family genes were found to be differentially regulated in this study. SLC8A1 and SLC43A2 were found to be up-regulated, whereas SLC39A10 was found to be down-regulated in the HUSFA group (Table 2). Two members of SLC genes (SLC16A7 and SLC27A6) were reported to be involved in FA metabolism [16]. Kaler and Prasad [31] postulated that SLC39A10 plays an essential role in cell proliferation and migration. However, the mechanism of SLC39A10 downregulation in FA metabolism is not yet clear, so further investigations are warranted to elucidate the function of this novel transcript regarding to FA metabolism. Sodhi et al. [32] reported that Glutamine fructose-6-phosphate transaminase 1 (GFPT1) is involved in glucose metabolism and differentially expressed in adipose tissue. A mutation in the exon of LEPR (p. Leu663Phe) is reported to be associated with increased feed intake and fatness in pigs [33].
Another gene family found to be differentially expressed that includes CYP17A, GSTCD and FABP7. These three genes were found to be down regulated in the higher USFA sheep in this study. Cytochrome P450 17A1 (CYP17A1, 17α-hydroxylase, 17,20-lyase) belongs to the cytochrome P450 super family that is expressed in the adrenals and gonads [34]. CYP2A6 gene is reported to be involved in meat flavour and odour-related molecules metabolism in sheep [35]. Barone et al. [36] reported that overexpression of CYP17A1 mRNA is associaed with enhancement of conjugated linoleic acid (CLA). The CLA refers to a group of positional and geometrical isomers of linoleic acid (cis-9, cis-12-octadecadienoic acid), an omega-6 essential fatty acid, that exhibit various physiological effects including anti-adipogenic, anti-carcinogenic, and immunomodulatory effect [37]. Glutathione S-transferase, C-terminal domain (GSTCD) belongs to the Glutathione S-transferases (GSTs) family that are functionally diverse enzymes, mostly known to catalyse FA conjugation reactions [38]. The GSTs transport different molecules [38] imply that GSTCD might transport FA to the tissues and thus involved in the FA metabolism in sheep. This study found that genes playing roles in fatty acid-binding protein (FABPs) were deregulated in higher USFA samples. Fatty acid-binding proteins such as B-FABP or FABP7 are known to be involved in the intracellular transport of PUSFA [39]. FABPs are intracellular proteins involved in binding and intracellular trafficking of FA for metabolism and energy production [40].

Biological function analysis for DEGs
Functional analysis showed that GO categories: biological processes, cellular components, and molecular functions were enriched in this study (Fig 3). The enriched biological processes identified were mainly related to cytokinesis, glycoprotein metabolic process, mitotic spindle, protein N-linked glycosylation, acute inflammatory response, and regulation of developmental process. Mitotic spindle organization plays a role in FA metabolism and energy productionin mammalian cells [41]. Cellular components consisted of cell projection part, extracellular space, integral to plasma membrane, and proteinaceous extracellular matrix were significantly enriched by the DEGs. Among the cellular components, proteinaceous extracellular matrix plays a role in skeletal muscle development in wagyu cattle [42]. The molecular functions identified were mostly related to kinase inhibitor activity, growth factor binding, GTPase activity, carbohydrate binding. It has been reported that growth factor binding is associated with serum insulin-like growth factor binding, thus influence lipid composition [43]. Carbohydrate binding is an important factor that influences FA metabolism in rat [44].
A total of 11 significantly enriched KEGG pathways were identified for DEGs (Fig 4). Pathway analysis revealed that glycosaminoglycans biosynthesis-keratan sulphate (KS), adipokine signaling, galactose metabolism, endocrine and other factor-regulated calcium metabolism, mineral metabolism, and PPAR signaling pathways have important regulatory roles in FA metabolism in the liver tissues. Keratan sulphate plays a crucial role in cells growth, proliferation, and adhesion [45]. Adipokine signaling acts as a bridge between nutrition and obesityrelated conditions [46]. Galactose metabolism is important for foetal and neonatal development as well as for adulthood [47]. Endocrine and other factor-regulated calcium metabolism, and mineral metabolism pathways are involved in intracellular mineral and calcium transportation, thus play roles in muscle muscle growth. Other important over-represented pathways in higher USFA group were phagosome and PPARs signaling pathway which were previously reported to be responsible for amino acid metabolism in cattle [16]. Several genes (APOA5, FABP7 and CPT1C) belonging to PPAR signaling pathway are identified in this study which could be involved in the FA metabolism in the seep. Berger and Moller [48] reported that PPARs are nuclear hormone receptors that are activated by FA and their derivatives, and regulate adipose tissue development and lipid metabolism in skeletal muscle. PPAR alpha is known to be involved in lipid metabolism in the liver and skeletal muscle, as well as in blood glucose uptake [49,50]. The PPAR signaling pathway was identified as the most significantly over-represented pathway involved in FA composition in cattle using RNA-seq [16], suggesting that PPAR could have a key role in controlling FA metabolism in sheep.

Regulatory hub genes of the hepatic transcriptome network
Regulatory hub genes of the hepatic transcriptome network identified several key genes including SOCS3, CBX6, MCM4, ITGB3, TGFBR2, GPRASP1, CELSR3, SDC3, SPOCK1, SEL1L and LEPR, which were upregulated in the liver tissues with higher USFA sheep ( Fig  5A). The SOCS3 negatively regulates JAK2/STAT5a signaling, thus inhibits FA synthesis in cow [51]. ITGB3 gene affects marbling development by promoting lipid accumulation and facilitates hepatic insulin [52]. The potential downregulated Hub genes identified were ACTA2, GPRASP1, TPM2, TGM3, PTK6, and LTF (Fig 5B). ACTA gene controls muscle filaments and energy utilisation in muscle [53]. GPRASP1 is involved in Calcium (Ca2+) release by skeletal muscle [54]. We, therefore, speculated that the potential network hubs identified in this study might play important roles in the FA composition in sheep. The co-expression network illustrated that RACGAP1, MCM4, SDC3, CKAP2, RNASE6, PREX1, QSOX1, and FUT11 were the upregulated Hub genes (Fig 6A). RACGAP1 gene involved in oxidative functions in skeletal muscle cells [55]. QSOX1 gene is reported to be involved in meat quality, lipid metabolism, and cell apoptosis, and suggested to use as a biomarker for cattle breeding for superior meat quality [56]. The co-expression network illustrated that NRN1, TPM2, CDC42EP5, SSC5D, GPRASP1, and HRC were the downregulated Hub genes (Fig 6B). NRN1 gene was expressed in various mammalian tissues including lipid rafts of cell membrane [57]. TPM2 gene is reported to be involved in muscle marbling development and suggested to be a candidate gene for meat quality traits in cattle [58]. Although, most of the co-expression networks were individually involved in FA composition traits, however, they exert functions through participating in different directions which implies that the FA composition is influenced by gene expression changes, and it is a complex physiological process.

Association between candidate markers and phenotypes
Selected polymorphisms within the APOA5, CFHR5, TFGBR2, and LEPR genes were found to be associated with the fatty acid composition phenotypes in this study ( Table 6). The APOA5 is mapped on the ovine chromosome 15, which is an important factor for triglyceride rich lipoprotein (TLR) regulation [59]. A member of APO gene family, APOV1 also known as APOVL-DLII, is found to be down regulated in higher (UFA) sheep. This gene was previously reported to be associated with UFA in chicken [60]. Significant association between the variants in APOA5 gene and high triglyceride levels and FA composition have been previously documented in sheep [61,62]. APOA5 is expressed in the liver, and controls VLDL binding (very low-density lipoprotein) to lipoprotein lipase (LPL) during FA synthesis in skeletal muscle and adipose tissue [63]. The CFHR5 is a 65 kDa plasma protein, binds with C3b, a C-reactive protein. Transforming growth factor beta receptor member familly 2 (TGBR2) is a member of the TGF-beta signaling pathway, which is involved in many cellular processes including cell growth, differentiation, and cellular homeostasis in animals [16]. The TGBPR2 gene is reported to be involved in myristoleic (C14: 1) FA metabolism [64]. Leptin receptor (LEPR) is an adipocytokine that regulates energy intake and uses in animals. Note, these polymorphisms are novel in sheep, and no association study with meat quality traits and FA compositions was conducted previously, so it is difficult to compare the results of this study with previous research. The LEPR was reported to be significantly associated with saturated FA, monounsaturated FA and polyunsaturated FA in pigs [1,65]. The upregulation of LEPR in higher polyunsaturated FA group and significant association indicate that this gene and marker may control the FA metabolism in sheep. Therefore, it could be postulated that LEPR, as a putative candidate gene plays crucial role in regulating fatty acid composition and metabolism in sheep.

Conclusion
The hepatic whole genome expression signature controlling unsaturated fatty acids (FA) levels in the sheep meat is, to our knowledge, deciphered for the first time. RNA-Seq provided a high-resolution map of transcriptional activities in the sheep liver tissue. The improvements in sheep genome annotations may lead to better coverage and detailed understanding of genomics controlling FA metabolism. This transcriptome analysis using RNA deep sequencing revealed potential candidate genes affecting FA composition and metabolism. This study suggested that candidate genes such as as APOA5, SLC25A30, GFPT1, LEPR, TGFBR2, FABP7, GSTCD, and CYP17A might be involved in the hepatic FA metabolism, thus control FA composition in muscle. Furthermore, number of SNPs were detected in the hepatic DEGs, and their associations with muscle FA compositions were validated. This transcriptome and polymorphism analyses using RNA Seq combined with association analysis showed potential candidate genes affecting FA composition and regulation in sheep. It is speculated that these polymorphisms could be used as markers for FA composition traits. However, further validation is required to confirm the effect of these genes and polymorphisms in other sheep populations.

Animals and phenotypes
Tissue samples and phenotypes were collected from the Indonesian Javanese thin-tailed sheep. All sheep (n = 100) were slaughtered in PT Pramana Pangan Utama, IPB University, and used for phenotyping as well as for association analysis. Animal's breeding, rearing and management, growth performance, carcass and meat quality data were collected according to guidelines of the Indonesian performance test. Animals were slaughtered with an average age of 12 months, and 30 kg of liveweight in slaughterhouse, in accordance with the Indonesian Inspection Service procedures and was approved by the 'Institutional Animal Care and Use Committee (IACUC)" issued by IPB University (approval ID: 117-2018 IPB). Tissue samples from the longissimus muscle (at least 500g between the 12/13th ribs) of each animal (left half of the carcass) were removed for this study. Tissue samples from the longisimuss muscle and the liver were collected, frozen in liquid nitrogen immediately after slaughter and stored at -80˚C until used for RNA extraction. Similar tissue samples were collected and stored at -20˚C for FA analysis. Fatty acids (FA) compositions were determined for each sample using the extraction method regularly performed in our Lab following Folch et al. [66]. Briefly, muscle samples (~100 g) were grinded for FA composition. The lipids were extracted by homogenizing the samples with a chloroform and methanol (2:1) solution. NaCl at 1.5% was added so that the lipids were isolated. The isolated lipids were methylated, and the methyl esters were prepared from the extracted lipids with BF3-methanol (Sigma-Aldrich, St. Louis, MO, USA) and separated on a HP-6890N gas chromatograph (Hewlett-Packard, Palo Alto, CA, USA) as described previously [67]. Gas-chromatography/mass spectrometry (GC-MS) method was applied for the quantification of FA compositions [66,67]. The average of USFA (MUSFA and PUSFA) and SFA value for these selected animals were 30.60 ± 10.12 and 39.73 ± 9.22 μg/g, respectively. Sheep having average USFA �45.59% μg/g and �25.84% μg/g were considered as higher-USFA (HUSFA) and lower-USFA (LUSFA) group, respectively (Table 1). In case of SFA, sheep having a SFA level �23.92% and �44.69% were considered as lower-and higher-SFA samples, respectively. However, for the transcriptome study, six sheep with divergently higher (n = 3) and lower (n = 3) USFA levels were selected from the total sheep (n = 100) population (Table 1). Total RNA was extracted from liver tissues using RNeasy Mini Kit according to the manufacturer's recommendations (Qiagen). Total RNA was treated using one-column RNase-Free DNase set (Promega), and quantified using a spectrophotometer (NanoDrop, ND8000, Thermo Scientific). RNA quality was assessed using an Agilent 2100 Bioanalyser and RNA Nano 6000 Labchip kit (Agilent Technologies).

Library construction and sequencing
RNA integrity was verified by Agilent 2100 Bioanalyser1 (Agilent, Santa Clara, CA, USA), where only samples with RIN > 7 were used for RNA deep sequencing. A total of 2 μg of RNA from each sample was used for library preparation according to the protocol described in Tru-Seq RNA Sample Preparation kit v2 guide (Illumina, San Diego, CA, USA). RNA deep sequencing technology was used to obtain the transcriptome expression. For this purpose, fulllength cDNA library was constructed from 1 μg of RNA using the SMART cDNA Library Construction Kit (Clontech, USA), according to the manufacturer's instructions. Libraries of amplified RNA for each sample were prepared following the Illumina mRNA-Seq protocol. The prepared libraries were sequenced in an Illumina HiSeq 2500 as single-reads to 100 bp using 1 lane per sample on the same flow-cell (first sequencing run) at Macrogen Inc, South Korea. The sequencing data have been deposited in NCBI (Accession: PRJNA764003, ID: 764003). All sequences are analysed using the CASAVA v1.7 (Illumina, USA).

Differential gene expression analysis
According to the FA concentration, animals were divided into two divergent phenotype value group (HUSFA and LUSFA) to identify differential expression genes (DEGs). The differential gene expression analysis was designed to contrast the differences in the expression of genes between two divergent sample group. The R package DESeq was used for the DEG analysis with raw count data [68]. The normalization procedure in DESeq handles the differences in the number of reads in each sample. For this purpose, DESeq first generates a fictitious reference sample with read counts defined as the geometric mean of all the samples. The read counts for each gene in each sample is divided by this geometric mean to obtain the normalized counts. To model the null distribution of computed data, DESeq follows an error model that uses a negative binomial distribution, with the variance and mean associated with regression. The method controls type-I error and provides good detection power [68]. After analysis using DESeq, DEGs were filtered based on p-adjusted value 0.05 and fold change � 1.5 [69]. Additionally, the gene expression data was analysed using a Generalized Linear Model (GLM) function implemented in DESeq to calculate both within and between group deviances. As sanity checking and filtration step, we cross-matched the results from both analysis (padjusted � 0.05 and fold change � 1.5 criteria, and GLM analysis) and only those genes which appeared to be significant in both of the tests (p value � 0.05) were selected for further analysis.

GO and pathways analysis
For biological interpretation of the DEGs, the GO and pathways enrichment analyses were performed using the NetworkAnlayst online tool [70]. For GO term enrichment, we used the GO database (http://geneontology.org/) and for pathways enrichment we used Kyoto Encyclopedia for Genes and Genomes (KEGG) database (https://www.genome.jp/kegg/pathway.html) incorporated in the NetworkAnlayst tool. The hypergeometric algorithm was applied for enrichment followed by Benjamini and Hochberg (H-B) [74] correction of multiple test.

Network enrichment analyses
To identify the regulatory genes, the sub-network enrichment analysis was performed using the NetworkAnlayst online tool [70]. The tissue-specific protein-protein interactions (PPI) data from DifferetialNet Basha et al. [71] databases incorporated with NetworkAnalyst with medium percentile were used for the creation of liver specific PPI network. The orthologous human symbol of the DEGs were uploaded into the NetworkAnalyst to construct the liver tissue-specific PPI network. The default network created one bigger subnetwork "continent", and 14 smaller subnetwork "islands". All the islands contain only single seed gene; therefore, those were not considered further. For high performance visualization, the continent subnetwork was modified by using the minimize function of the tool. The network was depicted as nodes (circles representing genes) connected by edges (lines representing direct molecular interactions). Two topological measures such as degree (number of connections to other nodes) and betweenness (number of shortest paths going through the node) centrality were taken into account for detecting highly interconnected genes (hubs) of the network. Nodes having higher degree and betweenness were considered as potentially important network hubs in the cellular signal trafficking. In addition, liver specific genes co-expression networks were also constructed using the TCSBN database Lee et al. [72] incorporated into NetworkAnalyst tool.

Quantitative Real Time PCR (qRT-PCR)
The cDNA was synthesised by reverse transcription PCR using 2 μg of total RNA, SuperScript II reverse transcriptase (Invitrogen) and oligo(dT)12 primer (Invitrogen). Gene specific primers for the qRT-PCR was designed by using the Primer3 software [73]. In each run, the 96-well microtiter plate was contained each cDNA sample, and no-template control. The qRT-PCR was conducted with the following program: 95˚C for 3 min, and 40 cycles: 95˚C for 15 s/60˚C for 45 s on the StepOne Plus qPCR system (Applied Biosystem). For each PCR reaction, 10 μl iTaqTM SYBR1 Green Supermix with Rox PCR core reagents (Bio-Rad), 2 μl of cDNA (50 ng/μl) and an optimized amount of primers were mixed with ddH 2 O to a final reaction volume of 20 μl per well. All samples were analysed twice (technical replication), and the geometric mean of the Ct values were further used for mRNA expression profiling. The housekeeping genes GAPDH and β-Actin were used for normalization of the target genes which were previously used for similar purpose in sheep tissues by our group [20]. The delta Ct (ΔCt) values was calculated as the difference between the target gene and geometric mean of the reference genes: (ΔCt = Ct target −Ct housekeeping genes) as described in Silver et al. [74]. The final results were reported as the fold change calculated from delta Ct-values.

Gene variation analysis
For gene variation analysis, SNP calls were performed on the mapping files generated by TopHat algorithm using 'samtools mpileup' command and associated algorithms [75]. Of the resulting variants, we selected the variants with a minimum Root Mean Square (RMS) mapping quality of 20 and a minimum read depth of 100 for further analyses. The selected variants were cross-checked against dbSNP database to identify mutations that had already studied. We also crosschecked and filtered the variants by the chromosomal positions of these variants against DEGs, and retained only those variants which mapped to DEG chromosome positions in order to find out the differentially expressed genes that also harboured sequence polymorphisms. By this way, we were able to isolate a handful of mutations that mapped to DEGs from many thousands of identified potential sequence polymorphisms. Furthermore, in order to understand whether these identified polymorphisms were segregated either in only one sample group (higher USFA and lower USFA) or in both groups (higher and lower USFA group), we calculated the read/coverage depth of these polymorphisms in all the samples [76]. The identified SNPs were classified as synonymous or non-synonymous using the GeneWise software (http://www.ebi.ac.uk/Tools/psa/genewise/ last accessed on 20.04.2021) by comparing between protein sequence and nucleotides incorporated SNP position [77].

Validation of SNP and association study
For the validation of association study, a SNP in each of four highly polymorphic DEGs (APOA5, CFHR5, TGFBR2 and LEPR) as well as the genes to be played key role in the fatty acid metabolism were selected for association study (Table 6). A total 100 sheep were slaughtered, and the blood sample were taken for DNA extraction until we got a final concentration of 50 ng/ml DNA. The genotyping process were performed by PCR-RFLP (Polymerase Chain Reaction-Restriction Fragment Length Polymorphism) method. The PCR were performed in a 15 ml volume containing 1 ml of genomic DNA, 0.4 μl of primers, 6.1 μl of MyTaq HS Red Mix, 7.5 μl of nuclease water. The PCR product was checked on 1.5% agarose gel (Fischer Scientific Ltd) and digested by using the appropriate restriction enzyme. Digested PCR-RFLP products were resolved in 2% agarose gels. Effect of genotypes on fatty acid composition was performed with PROC GLM using SAS 9.2 (SAS Institute Inc, Cary, USA). Least square mean values for the loci genotypes were compared by t-test, and p-values were adjusted by the Tukey-Kramer correction [78].
Supporting information S1 Table.